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We present an algorithm that prepares thermal Gibbs states of one dimensional quantum systems 
on a quantum computer without any memory overhead, and in a time significantly shorter than 
other known alternatives. Specifically, the time complexity is dominated by the quantity jv"^"/ T , 
where N is the size of the system, \\h\\ is a bound on the operator norm of the local terms of the 
Hamiltonian (coupling energy), and T is the temperature. Given other results on the complexity 
of thermalization, this overall scaling is likely optimal. For higher dimensions, our algorithm lowers 
the known scaling of the time complexity with the dimension of the system by one. 



Many open problems in condensed matter physics con- 
cern strongly correlated quantum many-body systems. 
These are typically not solvable analytically, and we 
have to resort to numerical simulations. Unfortunately, 
numerical methods tend to fail for general hamiltoni- 
ans on these systems, due to the exponential scaling of 
the dimension of the corresponding Hilbert space. This 
problem is one of the main motivations for the quest 
of quantum computers. Indeed, quantum computers 
can efficiently simulate unitary evolutions of quantum 
many-body systems with local interactions [TJ|2], because 
they can inherently deal with exponentially large Hilbert 
spaces. 

Nevertheless, the preparation of the desired initial 
state is still a difficult problem in general j3HZ]. There 
have been several proposals to tackle this problem [8]- 
[T2] , Some significant alternatives have worse complex- 
ity scaling than ours [5J [11] , while others apply to a re- 
stricted set of systems [lOj [12] . The quantum metropolis 
algorithm [9 , in particular, might often be faster, but 
lacks complexity bounds. The classical algorithm pro- 
posed in [13] can be used to prepare ID quantum thermal 
states with only a polynomial time complexity overhead 
with respect to our method, but its (classical) memory re- 
quirements scale exponentially with inverse temperature, 
and it does not extend to higher dimensional systems. 

The time complexity of our method for one dimen- 
sional systems is dominated by the quantity N^ h ^ T , 
where TV is the number of subsystems, T is the tempera- 
ture, and 1 1 ft 1 1 is a bound on the operator norm of the local 
terms of the Hamiltonian, the interaction strength. Note 
that this scaling is polynomial in N. The memory of the 
quantum computer scales simply with A/", an exponential 
improvement over general classical algorithms. Our algo- 
rithm can also be massively parallelized, and when run in 
a cellular automaton architecture the memory scales as 
7Vll^H/ T , but the total time would be linear in N (the to- 
tal number of steps would still be the same). In two and 
higher dimensions, our method lowers the number of ef- 
fective dimensions by one. This results in an exponential 
speedup, but the exponential scaling with N remains. 

The overall scaling appears to be optimal: the known 
complexity of thermalizing ID quantum systems makes a 



guaranteed polynomial scaling with temperature highly 
unlikely [6j[7]. We also expect the grouping of ||ft||/T in 
the exponent by dimensional analysis. In other words, 
the relevant temperature scale is set by the Hamiltonian. 

It is easier to introduce this method by explaining the 
proposal in [IT] first. In order to prepare a thermal state 
of a given Hamiltonian, the probability of each eigen- 
state needs to be set to the correct Gibbs probability. 
We can encode the correct probabilities as amplitudes of 
a marked state of an ancillary system. If a projective 
measurement of the ancilla returns the marked state, we 
succeeded in preparing the target Gibbs state. The suc- 
cess probability goes like the probability of projecting a 
random state into the target Gibbs state. If we use the 
totally mixed state as the initial state of these projec- 
tions, the success probability is Z/d N , where d is the 
dimension of each local subsystem. As a result, the num- 
ber of trials scales exponentially with the system size. It 
is possible to obtain a quadratic speedup over this scal- 
ing using Grover's amplitude amplification [11 , but the 
algorithm still scales exponentially with system size. 

We overcome the problem of exponential time cost by 
dividing the overall procedure into a sequence of projec- 
tions and arranging them so that we only need to rebuild 
a small section after most failures (see Fig. [I]). We first 
thermalize small regions, which we merge recursively un- 
til we have thermalized the whole system. Only when 
the failures are close to the end of this recursive proce- 
dure do we need to rebuild big sections. A careful error 
analysis shows that each of these merging operations can 
be implemented with a cost independent of the system 
size and the quantum correlation length, resulting in a 
running time that is only polynomial in the system size 
and independent of the correlation length. This method 
trivially generalizes to higher dimensions and reduces the 
scaling of the cost with the system dimension by one com- 
pared to a direct projection. 

We implement the merging perturbatively. Assume 
that we are given access to copies of p oc e~^ H (from 
previous steps). The Hamiltonian H corresponds to the 
halves to be merged, but the procedure is more general. 
We want to generate the state p^ oc e -P( H + h ) ^ where 
ft corresponds to the link between the two halves. We 
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FIG. 1. The procedure to thermalize an 8-qubit chain. After 
thermalizing individual qubits at level k = 0, we pair them 
up and merge them by adding the Hamiltonian that connects 
the two qubits. This procedure is then repeated recursively as 
we merge two already thermalized regions of size 2 k at level 
k to obtain a thermalized chain of size 2 k+1 at level k + 1. 
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FIG. 2. The conjugation circuit. First the quantum phase 
estimation circuit (QPE) writes the energy onto an ancilla 
register. Then the unitary U rotates the second ancilla con- 
ditioned on the energy. We measure the second ancilla and 
restart if we do not obtain |0). 



will see how to generate (with high probability) the state 



-P{H+eh) 



We then repeat the same process to 



produce the sequence 



p = p (0) 



p W _> p (2c) 



(1) 



Every transformation in the sequence has some probabil- 
ity of failure, in which case we restart. If all of the steps 
succeed, we approximate the state p^ with an error of 
O(e). It is important to remark that all the errors in this 
paper are in the trace norm. That is, in this case, for in- 
put p, we build a state a such that \\a — p^\\tt G O(e). 

We update the state p oc e~^ H to p (e) ex e -^ H+eh "> , to 
first order in e, in two steps. The first step is probabilistic 
and updates the probabilities of the Gibbs state through 
post-selection. If it fails, we will be forced to restart, 
and this will be the dominant part for the cost of the 
algorithm. In the second step we update the eigenbasis 
to the eigenbasis of p( € \ We now assume perfect phase 
estimation and perfect dephasing, and later account for 
the cost and errors of these operations. 

The probabilistic transformation of the first step is a 
conjugation with e -e/3 V 2 , to first order in e. We assume 
that h > 0. We use phase estimation and post-selection 
as in Fig. |2j This is similar to the procedure in [TT] . 
The phase estimation is of the unitary e 27rz/l£ , with 1/t > 



h>0. It implements the map ^2 a P a \tE a )(0\, where P a 
is a projection of the system onto the eigenspace of h 
with energy E a . This energy gets written in an ancilla 
register initialized to |0). We rotate a second ancilla to 
(1 — e/3E a /2)\0) + • • • |1), conditioned on the value of the 
previous ancilla register (unitary U in Fig. |2j), to get the 
ma P Ea( X - ef3E a /2)P a \tE a 0}(00\ + • • • . We then undo 
the phase estimation, obtaining the map 



5^(1 - 6/3^/2)^100X001 + - - ■ 

a 

= (1 - e/3/i/2)|00)(00| 



(2) 



Finally we measure the ancilla, and fail unless we obtain 
|0). We denote the result obtained by applying the above 
map to p as 

Pprob (x (1 - e(3h/2)p(l - e/3h/2) . (3) 
The success probability is 

p>l-ep\\h\\- (4) 

It can be seen that the conjugation with 1 — ef3h/2 of 
the previous paragraph updates all the probabilities of 
the eigenstates correctly to first order in e. It also imple- 
ments the appropriate transformation for the degenerate 
subspaces of p. Next, we update the eigenbasis from H 
to H + e/i, using the adiabatic approximation. This ap- 
proximation can be seen as a consequence of the Zeno 
effect, which can be achieved through measurements or 
dephasing [HI [15] . Here we implement the Zeno effect 
directly. More specifically, we use pure dephasing in the 
eigenbasis of H + eh (see [T5]). 

Denote by {P^} the projectors on the eigenbasis of 
H + eh. For input p, the operation J^ fce P^pP k e is pure 
dephasing on this eigenbasis. To qualify the effect of 
dephasing on p pro b> we start by writing the result of this 
ideal operation as [T5] 

yZ p k e PprobPk* = lim / dtg a (t)U e (t)p proh U^(t) , (5) 

/c e J 

where g a is the density function of a Gaussian dis- 
tribution with mean and standard deviation a, and 
U e (t) = e -*(^+ e/l ) t . This equality is easy to check di- 
rectly. We now use a Dyson series to second order for 
the time evolution: 



U e (t) = U (t) - ieA(t) - e 2 B(t) 
U (t) = e~ mt 



A(t) = f 
Jo 



dtt U (t - ii)M/ (ii) 



(6) 
(7) 

(8) 



B(t)= [ dti [ dt 2 Uo(t-t 1 )hUo(t 1 -t 2 )hU e (t 2 ) ■ 
Jo Jo 
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We then insert terms to first order in e from Eq. (|6| 
into Eq. The second (and higher, after expanding 

further) order terms can all be bounded by 0(e 2 ||/i|| 2 ) in 
the trace-norm. 

We want to compare the state obtained after the de- 
phasing to p^ , to first order in e. With this goal, we use 
the Dyson series in imaginary time to expand p^ . In this 
expansion we need intermediate states p{fi) = e~^ H / 'Z, 
where the partition function is always at inverse temper- 
ature [), Z = Tre - ^. We obtain, to second order, 



P 



p [e) =p-e / dp 1 p{p-p 1 )hp{p 1 ) (9) 
Jo 

+e 2 / dp, / dfapiP-PJhpiPi-fohpifa). 
Jo Jo 

The second and higher order corrections can be bounded 
by (9(e 2 /3 2 ||/i|| 2 ) in trace-norm. We assume that P > 1, 
so this error dominates the error G(e 2 \\h\\ 2 ) above. 

We aim to obtain the first order expression in Eq. ([9| 
with our procedure. The spectral decomposition of p 
is ^ k PkPki where P k is the projector to the subspace 
spanned by eigenstates with eigenvalue E k of H, and 
Pk = e~^ Ek j ' Z. The term proportional to e in the expan- 
sion of p( e ) in Eq. ([9| can be rewritten in terms of these 
projectors, up to normalization, as: 



| p U1D , ^ PihPk±PkhPi , 



l^k 



Ei — E k 



Finally, using H = E k Pk in Eq. ([5]), with the expan- 
sion (p|), gives the first order correction Eq. (10). 



Until now, we were considering the behavior of our al- 
gorithm assuming perfect phase estimation and dephas- 
ing, which is not realistic. We now account for the effects 
of the errors inherent in the two parts of our perturbative 
Hamiltonian update algorithm when implemented using 
only dynamical evolution with the Hamiltonian H + kh 
for < k < 1 . We quantify the total cost in terms of 
the evolution time. 

For the conjugation circuit, we use high precision phase 
estimation [TTJ [T6HT8] . The cost (evolution time with h) , 
for accuracy S and error e, scales as 0(log(l/e)/5). This 
implements the transformation 



E p «® E^i^+^iu pi 



(ii) 



with t\E a — Ef\ < 6 and q a < e. The t is the evolution 
time of the unitary e 27Tlht used for the phase estimation, 
with 1/t > h > as above. The term |£ a ) groups all the 
states of the ancilla register with a reading outside the 
goal accuracy 5. 

First consider the effect of the error term q a Pa\€a)- 
All the manipulations conserve the projectors P a and do 



not increase the norm of |£ a ). Therefore, the final error 
due to these terms, on input a, can be bounded with 
terms like || J2 a Qa,Pa&\\tT < £• The final error can in- 
crease as a result of the normalization when projecting 
onto the post-selected state, if the preparation is suc- 
cessful. We will see that this effect is negligible. We then 
rotate a second ancilla and undo the phase estimation to 
obtain an approximate version of the conjugation unitary 
in Eq. ([2| . The undo of the phase estimation incurs in an 
error bounded by 0(ep5) in the trace norm, and the rest 
of the errors are O(e) as before. If we choose e and ePS 
to be 0(e 2 p 2 \\h\\ 2 ), we can implement the map in Eq. (|3| 
with the same probability Q as before, trace norm error 
G(e 2 p 2 \\h\\ 2 ), and cost (evolution time) 



0(log(l/(e^||/i||))/(e/?||/ l || 2 )). 



(12) 



Finally, we deal with the errors related to imperfect de- 
phasing. Irrespective of the implementation of dephas- 
ing, imperfect accuracy in the dephasing results in an 
operation where only phases between eigenstates with 
relative gap bigger than some bound £ are suppressed. 
This is modeled by the operation 



E 

:|^-£jl<C 



PheaPr 



(13) 



(cf. Eq. (|5|). We obtain imperfect dephasing if we choose 
the standard deviation of the Gaussian in Eq. ([5| to be 
a = (9(1/C), instead of (approaching) infinity [19]. 

We analyze the effect of imperfect dephasing errors by 
the following purely formal mathematical procedure. For 
a given Hamiltonian H + e/i, we define the Hamiltonian 
H which is similar to H + eh, but such that all the gaps 
are at least £. We do this by increasing the degener- 
acy: we group all eigenvalues in bins with relative gap 
£. The imperfect dephasing (or accuracy () would be 
fundamentally exact if the target Hamiltonian was i7, 
because all the gaps would be bigger than the accuracy. 
Let x = H — (H + e/i), with ||x|| < C m operator norm. 

By an expansion similar to that of Eq. ([9| we can write 
p (e) = ~( e ) + O(PC), where p^ = e~P 6 /Tr and the 

bound is, as always, in the trace norm [20] . Similarly, we 
can write U e (t) in terms of U e (t) = e~ %Ht and x using 
Dyson series. Plugging these expressions into Eq. ([5| 
with a = c/£ for a constant c = 0(log(l/e)), we get 



/ 



dtg c/<: (t)U e (t)p pioh uHt) = p (£) + 0{fiQ , (14) 



where g c j^ is a Gaussian with standard deviation c/(. 



Now choosing £ = e 2 /3||/i|| 2 we see that we can transform 
p to p( e ) with the probability of Eq. Q , trace norm error 
G(e 2 p 2 \\h\\ 2 ), and cost (evolution time) 



0(log(l/(eP\\h\\))/(e 2 p\\hf)) 



(15) 
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We can merge two regions already thermalized into 
one large thermal region with the two subroutines just 
described using a sequence of small perturbative steps 
(see Eq. Q). Each step is successful with the proba- 
bility given in Eq. @. The average number of steps 
until we generate a complete sequence without failures 
is (m) G 0(e^) [21]. Each time that we fail we need 
to produce two new thermal regions to be merged. The 
average number of failures is (a) £ 0(e^ h ^). 

Now we analyze the average number of steps (r(k)) 
required to prepare a thermalized chain of length 2 k at 
level k (see Fig. [I]). Since a and r(k — 1) are independent 
random variables, the expectation value of r(k) is 



<r(fc)>=2(a><r(fc-l)> + <m>. 



(16) 



This gives (r(logiV)) £ <9(exp((/?||h||+log2)logiV)). Sim- 
ilarly, the total error is 0(Ne(3 2 \\h\\ 2 ). If we choose e = 
e/ (N fi 2 \\h\\ 2 ), we get a total error of 0(e) in trace-norm. 



Finally, using Eq. ( 15 ) for the evolution time of each step, 



we obtain the dominant contribution to the total evolu- 
tion time /3N^ h W/e 2 . 

We have presented an algorithm that prepares a ther- 
mal state of a ID quantum system in time polynomial 
in the system size and exponential in the inverse tem- 
perature (as required by the existence of Q MA- complete 
ground state problems in ID). This algorithm can be gen- 
eralized into D dimensions. At level k of the recursion, 
we have built squares (for 2D) or cubes that are now 
merged. We do not get polynomial scaling with system 
size for D > 1 because the intersection of two neigh- 
boring regions goes like TV^ -1 . Note that this is to be 
expected because there exist 2D ground states with con- 
stant gap that encode the solution to NP-complete prob- 
lems. A careful analysis confirms that the time complex- 
ity is dominated by the operations at the top level, and 
the dominating factor is /3e 2 ^ h ^ DN /e 2 . This is an 
exponential speedup from the known ex.p(0(N D )). The 
memory requeriments still scale with the number of sites 
of the model, N D . 

There are also several possible improvements to the 
scaling of this algorithm. If one is interested in thermaliz- 
ing a classical system with a small quantum perturbation 
one can first solve for the classical part of the Hamilto- 
nian. Then, one would only need to use projections for 
the quantum perturbation. Also, if one is interested in 
thermalizing a quantum system with short-ranged quan- 
tum correlations, one can also use belief propagation [22]- 
25J to reduce the storage requirements from O(N) qubits 
to 0(1 log(iV)), where I is a constant related to the quan- 
tum correlation length. This can be done by tracing out 
parts of the blocks that do not share any entanglement 
with the boundary to be merged. The time complexity 
of the algorithm remains the same as before. Note that 
the cost (both in memory and time) of the classical al- 



gorithm of quantum belief propagation for ID quantum 
systems is exponential in I ~ 1/T. This bound is only 
heuristic, and similar to [13]. 
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